# Code to create Figure 1 in "Pandemics and Cities: Evidence from the Black Death and the Long-Run"
# Author: Noel Johnson
# Last updated: 9-14-23

library(sf)
library(raster)
library(tidyverse)
library(tmap)

# set the working directory
setwd("/Users/noeljohnson/Dropbox/Research/Plague City Growth/JUE RR1 Nov 2022/figures_replication/Figure_1")

# load in city data
# define the projection as europe equidistant conic (EEC)
EEC <- "+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=43 +lat_2=62 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs"
# read in the city data
cities_165 <- read.csv("cities_for_map165.csv")
# spatialize the city data
cities_165_spatial <- st_as_sf(cities_165, coords = c("longitude", "latitude"), crs = 4326)
# project the cities data into EEC
cities_165_spatial_proj <- st_transform(cities_165_spatial, EEC)
# plot the cities to make sure they look correct
plot(cities_165_spatial_proj["mortality"], reset = FALSE)

# Load in State Borders
# read in the modern country borders shape file
modern_borders <- st_read("ModernCountries/Modern Europe Projected.shp")
# project the modern border into EEC
modern_borders_proj <- st_transform(modern_borders, EEC)
# look at the extent of the projected shape file
modern_borders_proj
# original Extent of modern countries c(xmin=-2639987, xmax=1655719, ymin=3243609, ymax=8097861))
# clip the modern countries shape file
modern_borders_clip <- st_crop(modern_borders_proj, c(xmin=-1100000, xmax=1655719, ymin=3243609, ymax=6800000))
# drop poland and slovakia from the map since we don't have any cities there
modern_borders_clip <- modern_borders_clip %>%
  filter(CntryName != "Poland") %>%
  filter(CntryName != "Slovakia")

# plot the border with the cities to make sure it looks correct 
plot(modern_borders_clip[6], reset = TRUE)
plot(cities_165_spatial_proj["mortality"], add = TRUE)

# make figure 1 in greyscale

# first make the modern country borders map...
state_map <- tm_shape(modern_borders_clip) + tm_borders(col = "black") 

# Define the breaks for mortality for both size and color
breaks <- c(25, 50, 75, 100)

# Custom grayscale color palette based on breaks
color_palette <- gray.colors(length(breaks) - 1) # because you have one less interval than breaks

# Extract coordinates from the cities sf object
cities_165_spatial_proj$lon <- sf::st_coordinates(cities_165_spatial_proj)[, 1]
cities_165_spatial_proj$lat <- sf::st_coordinates(cities_165_spatial_proj)[, 2]

# Construct the plot
combined_plot <- ggplot() +
  geom_sf(data = modern_borders_clip, color = "black", fill = NA) +
  geom_point(data = cities_165_spatial_proj, 
             aes(x = lon, 
                 y = lat,
                 size = mortality, 
                 fill = mortality),
             shape = 21,  # This shape supports fill
             color = "white",
             show.legend = TRUE) + 
  scale_size_continuous(name = "Mortality Rate",
                        breaks = breaks,
                        range = c(1, 6)) +
  scale_fill_gradient(name = "Mortality Rate",
                      low = "lightgrey",
                      high = "black") +
  guides(size = "none") +  # This line removes the size legend
  theme_minimal() +
  theme(
    legend.position = c(0.08, 0.45),     # Position the legend slightly to the right while keeping it centered vertically
    legend.key.size = unit(0.35, 'cm'), # Adjust the size of the legend keys
    legend.text = element_text(size = 10),  # Adjust the size of legend text
    legend.title = element_text(size = 11), # Adjust the size of the legend title
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),        # Remove axis labels
    axis.title = element_blank()        # Remove axis titles
  )

print(combined_plot)

ggsave(filename = "figure_1.pdf", plot = combined_plot, width = 10, height = 7, units = "in")



#




#End of Code



# # a map in color (not for the journal)
# # first make the modern country borders map...
# state_map <- tm_shape(modern_borders_clip) + tm_borders(col = "black") 
# state_map
# # now add the cities with color and size of circles representing mortality rate
# mycols <- colors()[c(142, 54, 35, 41, 24)]
# cities_map_color <- tm_shape(cities_165_spatial_proj) +
#   tm_bubbles("mortality", col = "mortality", border.col= "black", border.alpha = 0.5, style="fixed", breaks = c(-Inf, 20, 40, 60, 80, Inf), palette=mycols, title.size = "Mortality Rate (Size)", scale = 2.0, title.col ="Mortality Rate (Color)") +
#   tm_layout(legend.stack = "vertical", legend.text.size = .80, frame = F, inner.margins = .10, legend.position = c(-.045, 0.25)) +
#   state_map
# # view the map
# cities_map_color
# # save the map
# tmap_save(cities_map_color, filename="cities_map_color.png",
#           height=8.5, width=11, units="in", dpi=300)










